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ABSTRACT 

We present a new backward evolution model for galaxies and AGNs in the infrared 
(IR) . What is new in this model is the separate study of the evolutionary properties 
of the different IR populations (i.e. spiral galaxies, starburst galaxies, low-luminosity 
AGNs, "unobscured" type 1 AGNs and "obscured" type 2 AGNs) defined through a 
detailed analysis of the spectral energy distributions (SEDs) of large samples of IR 
selected sources. The evolutionary parameters have been constrained by means of all 
the available observables from surveys in the mid- and far-IR (source counts, redshift 
and luminosity distributions, luminosity functions). By decomposing the SEDs repre- 
sentative of the three AGN classes into three distinct components (a stellar component 
emitting most of its power in the optical/near-IR, an AGN component due to hot dust 
heated by the central black hole peaking in the mid-IR, and a starburst component 
dominating the far-IR spectrum) we have disentangled the AGN contribution to the 
monochromatic and total IR luminosity emitted by the different populations consid- 
ered in our model from that due to star-formation activity. We have then obtained 
an estimate of the total IR luminosity density (and star-formation density - SFD - 
produced by IR galaxies) and the first ever estimate of the black hole mass accretion 
density (BHAR) from the IR. The derived evolution of the BHAR is in agreement with 
estimates from X-rays, though the BHAR values we derive from IR are slightly higher 
than the X-ray ones. Finally, we have simulated source counts, redshift distributions 
and SFD and BHAR that we expect to obtain with the future cosmological Surveys 
in the mid-/far-IR that wiU be performed with JWST-MIRI and SPICA-SAFARI. 

Key words: cosmology: observations - galaxies: active - galaxies: evolution - galax- 
ies: Seyfert - galaxies: starburst - infrared: galaxies. 



1 INTRODUCTION 

In the past years strong observational evidence of high rates 
of evolution for infrared (IR) galaxies has been obtained 
by means of two independent findings: the detection of a 
large amount of energy contained in the Cosmic Infrared 
Background (CIRB, Hauser & Dwek 2001), and the num- 
ber counts from several deep cosmological surveys (from 15 
^m to 850 /im) largely exceeding the no-evolution expecta- 
tions. Both results agree in requiring a strong increase in the 
IR energy density between the present time and The 
discovery of the CIRB, which is interpreted as the integrated 
emission from dust present in galaxies, has offered new per- 
spectives on our understanding of galaxy formation and evo- 
lution, since it provides a key constraint on the history of 
star-formation (SF) and accretion in the Universe. The res- 
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olution of the CIRB into individual sources has been one 
of the main goals for the Infrared Space Observatory {ISO; 
Kessler et al. 1996) and the Spitzer Space Telescope (Werner 
et al. 2004). The deep cosmological surveys carried out by 
ISO and Spitzer allowed new insights into the IR population 
contributing to the CIRB, showing source counts in excess 
with respect to no-evolution predictions and revealing the 
existence of a population of distant, dusty IR-bright galaxies 
that are missed by optical surveys. The major contributors 
to the CIRB are Luminous (Lg-iooOfim > 10" Lq, LIRGs) 
and Ultra-Luminous IR Galaxies (Lg_iooonm > 10^^ Lq, 
ULIRGs): they have been discovered by the Infrared Astro- 
nomical Satellite {IRAS; Neugebauer 1984) as rare objects in 
the local Universe radiating most of their energy in the IR, 
but have been successively found by ISO, Spitzer and also 
SCUBA on the JCMT in the sub-millimeter (sub- mm), to 
become more and more important as the redshift increases. 
In particular, the so-called "sub-mm galaxies" (SMGs) de- 
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tected by SCUBA (i.e. Smail, Ivison & Blain 1997) and sub- 
sequently confirmed by Spitzer (i.e. Pope et al. 2006), which 
emit a significant fraction of their rest-frame bolometric lu- 
minosity in the far-infrared (FIR)/sub-mm and seem to be 
mostly at 2 ~ 2 — 3, represent a direct evidence for that. 
The existence of this strongly evolving population of dusty, 
massive galaxies that form the bulk of their stars at high 
redshifts (SMGs and high-z (U)LIRGs are found to be al- 
ready massive galaxies at 2 ~ 2; i.e. Dye et al. 2008), is an 
example of anti-hierarchical behaviour, showing how crucial 
this IR-bright, dust-obscured galaxy population is to under- 
stand galaxy formation and evolution. 

The new results provided by the Balloon-borne Large 
Aperture Sub-millimeter Telescope {BLAST; Pascale et al. 
2008) and the Herschel Space Observatory (Pilbratt et al. 
2010) in the FIR/sub-mm domain (e.g., Patanchon et al. 
2009; Bethermin et al. 2010; Berta et al. 2010; Oliver et al. 
2010; Gruppioni et al. 2010), together with the availability 
of the new space facilities in the coming years, such as James 
Webb Space Telescope {JWST; Gardner et al. 2006) and fur- 
therly the SPace Infrared telescope for Cosmology and Astro- 
physics {SPICA; Nakagawa et al. 2009), open a new perspec- 
tive to study in detail the population of IR galaxies beyond 
z = 1 — 2, requiring new models to explain the high rate of 
evolution observed at lower redshifts. Two main weaknesses 
exist in most of the current models for IR sources: 1) the 
failure in reproducing the observed redshift distributions; 
2) the severe underestimate of the contribution from Active 
Galactic Nuclei (AGNs). The new models should take into 
account that IR galaxies can host both SF and AGN activ- 
ity (i.e. Lutz et al. 1998) and that in the far-infrared (FIR) 
even the emission from Seyfert and Quasars can be domi- 
nated by SF (Lutz et al. 2004; Schweitzer et al. 2006; Shao 
et al. 2010). We can no longer neglect the AGN contribu- 
tion in modelling and interpreting IR data, but we should 
rather understand and quantify the AGN presence within IR 
galaxies and its connection with SF activity. In particular, 
a key cosmological question that needs to be answered by 
models and future observations regards the role of AGN in 
galaxy formation and evolution. In fact, the seminal discov- 
ery that all massive galaxies in the local Universe harbour 
super-massive black holes (SMBH; Mbh > IO'^Mq) imphes 
that all massive galaxies have hosted AGN at some time 
during their life (i.e. Magorrian et al. 1998). Recent works 
in the X-rays, optical and mid-infrared (MIR) have shown 
that many heavily obscured AGNs escape detection even in 
the deepest optical and X-ray observations, but, as expected 
from their high level of obscuration, reveal themselves in the 
IR (i.e. Daddi et al. 2007; Fiore et al. 2008). When planning 
future IR surveys with new facilities and/or when trying to 
interpret new data results, it is therefore necessary to con- 
sider models that properly take into account the presence 
of AGNs within a significant fraction of the IR population, 
in order to be able to answer the still open questions about 
galaxy formation and evolution, like i.e. What drives the 
evolution of the massive, dusty distant galaxy population? 
What feedback/interplay exists between the AGN and star- 
forming phase and how does it relate to cosmic downsizing? 

We have developed a new backward evolution model 
fitting all the main constraints provided by the IR/sub-mm 
surveys (from 15 fim to 500 t-tm), where we have taken 
particular care in properly identifying and modelling the 
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AGN contribution at diflerent luminosities and redshifts. 
Our model starts from the classical approach consisting of 
evolving a local luminosity function (LLF) in luminosity 
and/or density with redshift. In particular, we have decom- 
posed the LLF into diflerent populations of galaxies and 
AGNs, each of them evolving independently. What is re- 
ally new in this model is the way we have distinguished 
between the different IR populations: the separation into 
classes is based on a detailed Spectral Energy Distribution 
(SED) study performed on a large spectroscopic sample of 
MIR selected sources (Gruppioni et al. 2008). The shape 
of the SED provides a more meaningful physical character- 
isation of different populations than e.g. their luminosity, 
which is commonly used to separate IR sources into differ- 
ent classes (i.e. LIGs, ULIGs, HyLIGs). We have constrained 
our model by means of the MIR (i.e. from ISO and Spitzer) 
and FIR (i.e. from Herschel) data. In particular, as con- 
straints in the FIR we have considered the very recent re- 
sults of the PACS Evolutionary Probe (PEP) Survey (Berta 
et al. 2010; Gruppioni et al. 2010), of the _ffersc/ie/-ATLAS 
(H-ATLAS) Survey (Eales et al. 2010; Clements et al. 2010) 
and of the Herschel Multi-tiered Extra-galactic (HerMES) 
Survey (Oliver et al. 2010; Vaccari et al. 2010). The model 
evolutions and luminosity functions have then been used to 
derive the evolution with redshift of the star-formation and 
accretion densities from IR luminosity. To this purpose, we 
have analysed the relative contributions of AGN and star- 
formation activity to the IR luminosity of the populations 
considered in the model by means of a SED decomposition 
obtained with the Fritz et al. (2006) code. 

The paper is structured as follows. In Section[2]we pro- 
vide a detailed description of our model, including the pop- 
ulation separation, the considered MIR and FIR data con- 
straints, the LLF parameters and the evolutionary paths. 
Section |3] reports on the results of our model in the MIR 
and FIR/sub-mm bands, while in Section [4] we discuss the 
contribution of AGNs to the total IR luminosity and its evo- 
lution, deriving also an estimate of the Black Hole accretion 
rate density (BHAR) and star-formation density (SFD) evo- 
lution with redshift. In Section [5] the model predictions for 
Surveys with the Mid-InfraRed Instrument {MIRI; Wright 
et al. 2004) on board of JWST and with the SPICA FAR- 
Infrared Instrument (SAFARI; Swinyard et al. 2009) are dis- 
cussed, while our conclusions are presented in Section |6] 



2 THE MODEL 

The model starts from the classical approach of evolving 
a local luminosity function (LLF) with redshift, but con- 
siders five different IR populations evolving independently. 
As starting LLF we consider a parametrisation of the 12- 
^m one derived by Fang et al. (1998), consistently with the 
Rush, Malkan & Spinoglio (RMS, 1993) LLF, decomposed 
into single population LLFs as to reproduce the total ob- 
served LLF. As parametrisation for the LLF ($(!/), defined 
as the differential LF per decade in luminosity), we adopt 
the form first derived by Saunders et al. (1990) to describe 
the 60-/im LF of IRAS galaxies, found to be too broad to 
be described by a standard Schechter function. Saunders et 
al. (1990) found a better representation for the IR LF with 
"1>(I/) given by the function: 
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Figure 1. Model LLF: spiral galaxy (green dashed), starburst 
galaxy (cyan dot-dashed), LLAGN (red dot-dot-dot dashed), AGN2 
(magenta dotted), AGNl (blue long-dashed), total (black solid). 
The different populations are described in Section |2.1| The ob- 
served Fang et al. (1998) and RMS data points are overplotted 
as black filled circles and open squares respectively. 
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which behaves as a power law for L<<L* and as a Gaussian 
in logL for L>>L*. 

The relative fractions of sources within the different 
classes are defined on the basis of a detailed broad-band 
SED-fitting analysis performed on a large sample of MIR 
(15-^m) selected galaxies, all with spectroscopic redshift and 
classification (Gruppioni et al. 2008). The MIR sample on 
which we have originally based the source classification con- 
sists of 203 sources from the ELAIS-Sl survey at 15 /im 
(Lari et al. 2001; Gruppioni et al. 2002), for which a de- 
tailed optical spectroscopical analysis has been performed 
and presented by La Franca et al. (2004). Based on the SED- 
fitting technique Gruppioni et al. (2008) have classified the 
MIR sources, identifying AGN signatures in about 50% of 
them. Similar AGN fractions in IR surveys have been re- 
cently found either in local samples (i.e. Smith et al. 2008; 
Goulding and Alexander 2009) or at cosmological distances 
through different identification techniques (i.e. Brand et al. 
2006; Treister et al. 2006). The parametrisation of the single 
LLPs, as well as the evolution parameters for each popula- 
tion have been obtained through a minimisation algorithm 
aimed at finding the best-fitting parameters to simultane- 
ously reproduce the redshift and luminosity distributions 
observed for the five populations of the Gruppioni et al. 
(2008) sample (with the same optical and IR completeness 
corrections and redshift limits of the sample data applied 
to the model), the 15-fim to 500-/im source counts (as de- 
scribed in the following Sections). 

In Figure [l] we show the LLF decomposition into the 
five different populations (whose parameters are reported in 
Table [l]), overplotted to the Fang et al. (1998) and RMS 
12-//m observed data points. 

2.1 The Infrared Populations 

As mentioned in the previous section, the classification of the 
IR sources into different populations in our model is based 
on the results of a broad-band SED-fitting analysis for one of 
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Figure 2. Example of an observed SED (red filled circles) of a 
MIR, source spectroscopically classified as galaxy, compared to 
three different template SEDs (from Polletta et al. 2007): Seyfert 
2 (black solid line, best-fit), Starburst Galaxy (green dashed) and 
Sd Spiral Galaxy (blue dot-dashed), normalised to the optical 
band. 



Table 1. 15-/^m LLF Parameters 



SED class 


logio(L*(0)) 


<j>* 


Q 




spiral 


9.45 


1.7x10^3 


1.35 


0.30 


starburst 


9.70 


5.0x10-^ 


0.05 


0.31 


LLAGN 


9.76 


4.3x10-* 


0.99 


0.18 


AGN2 


9.95 


3.0x10"^ 


0.001 


0.27 


AGNl 


9.80 


4.0x10-5 


1.65 


0.60 



the largest available spectroscopic samples of MIR-selected 
galaxies and AGNs at intermediate redshifts (z < 1.5) per- 
formed by Gruppioni et al. (2008). The sample, consisting 
of 72% of the 15-Atm ELAIS-SWIRE Survey in SI (Lari et 
al. 2001; Gruppioni et al. 2002), contains 203 extragalactic 
sources, all with measured spectroscopic redshift (La Franca 
et al. 2004). The sample of 203 15-/im sources is composed 
by a deeper sample in the central Sl_5 area, which is 99% 
spectroscopically complete to R — 21.6, plus a shallower 
sample in the rest of the field, which is 97% spectroscopi- 
cally complete to _R = 20.5. The same optical/15-pim limits 
of the data have been applied to the model when comparing 
the model to the data in ELAIS-Sl (see the single popu- 
lations source counts, redshift and luminosity distributions 
plotted in Figures 5 and 6). Most of these sources have full 
multi- wavelength coverage from the far-UV (GALEX; Mar- 
tin et al. 2005) to the FIR (Spitzer) and lie in the redshift 
range 0.1 < z < 1.3. This large sample allowed us for the 
first time to characterise the spectral properties of sources 
responsible for the strong evolution observed in the MIR 
and to construct an observational library of templates for IR 
galaxies and AGNs at intermediate z. The observed SEDs 
have been interpreted by performing a fit with several lo- 
cal template SEDs representative of different classes of IR 
galaxies and AGNs (Polletta et al. 2007), comparing the re- 
sulting SED classification with the spectroscopic one. The 
considered templates included three elliptical galaxies of dif- 
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Figure 3. Template SEDs considered for the five IE, populations. Top left: spiral SEDs, evolving with Lis^m from an Sa spiral to an 
Sdm one, with fine interpolations in between. Top right: starburst SEDs evolving from the moderate starburst M82 to the extreme 
stajrburst Arp220. Middle left: Seyfert 2 SEDs (original PoUetta's and modified - e.g. with a higher FIR bump - as to reproduce a 
fra<;tion of Herschel sources showing an excess in the FIR with respect to the original template), assumed as representative of the LLAGN 
population. Middle right: Markarian 231 SED, assumed as template for the AGN2 population. Bottom right: TQSOl SED, assumed for 
the AGNl population. 



ferent ages, one lenticular, seven spirals, three starbursts, 
three type 1 QSOs, one type 2 QSO, Seyfert 1, 1.8 and 2 
and two composite ULIRGs, containing both starburst and 
AGN component, in the wavelength range between 0.1 and 
1000 /im. Of the original sample of 203 somxcs 41% is well 
reproduced by galaxy templates (SO to Sdm, no ellipticals), 
6.5% by starburst templates, 35% by Seyfert 2/Seyfert 1.8 
templates, 5.5% by templates typical of "obscured" AGNs, 
characterised by large column densities of obscuring mate- 
rial (i.e. Markarian 231 or IRAS19254) and 12% by "unob- 
scured" AGN templates. Note that the two SEDs reproduc- 
ing those of the the "obscured" AGN population are empiri- 
cal templates created to fit the SEDs of the heavily obscured 
broad absorption-line (BAL) QSO Markarian 231 (Berta 



2005) and the Seyfert 2 galaxy IRAS 19254-7245 South 
(Berta et al. 2003). These two SEDs are similar in shape, 
containing a powerful starburst component, mainly respon- 
sible for their FIR emission, and an AGN component that 
contributes to the MIR (Farrah et al. 2003), and reproduce 
the SEDs of "obscured" AGNs regardless of their optical 
spectra (i.e. broad or narrow lines in the optical; Gruppioni 
et al. 2008). The fraction of objects showing different levels 
of AGN activity in their SEDs is significantly higher (~53%) 
than that derived from the optical spectroscopy (~29%; La 
Franca et al. 2004), particularly because of the identifica- 
tion of AGN activity in the SEDs of objects spectroscopi- 
cally classified as galaxies. This might be partially due to 
the fact that the spectroscopic classification can be some- 
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what unreliable because of dilution of AGN signatures due 
to the host galaxy light in the optical band. It is likely that 
in most of these objects the AGN is either obscured or of 
low-luminosity, and thus it does not dominate the energetic 
output at any wavelength, except in the MIR, showing up 
just in the range where the host galaxy SED has a minimum. 
Similar results have been obtained from the same detailed 
broad-band SED-fitting analysis performed on a sample of 
100- and 160-^m Herschel sources in the PEP GOODS-N 
field (Gruppioni et al. 2010). In the fields where we have 
performed detailed SED-fitting analyses, we have done fur- 
ther investigations on the available X-ray images to check 
whether a correlation between our SED classification and 
the X-ray detection and luminosity exists or not. Indeed, 
we find that the X-ray detected sources are mainly objects 
classified as AGNs on the SED-fitting basis, with their lumi- 
nosities being higher than those of the few detected galaxy- 
SED ones. In particular, from a match between the posi- 
tions of IR selected sources classified by Gruppioni et al. 
(2008; 2010), with those of the X-ray sources in ELAIS- 
Sl (XMM: Puccetti et al. 2005; Chandra: Vignah et al. in 
preparation) and in GOODS-N (Chandra: Alexander et al. 
2005) respectively, we found that, indeed, the fraction of 
matches for the Seyfert2/1.8 SED-classified objects is sig- 
nificantly higher than those found for objects classified as 
starburst and spiral galaxies on the basis of their SEDs. In 
particular, in ELAIS-Sl we found 41% of XMM+Chandra 
detections for the Seyfert2/Seyfertl.8 sources against 15% 
for the starburst -|-spiral galaxies, while in the GOODS-N we 
found 25% of Chandra detections for the Seyfert 2's (with 
<L(2-8keV)>~10"*^ erg s"^) against 12% and 4% for the 
starburst (<L(2-8keV)>~10'*^'^ erg s~^) and spiral galax- 
ies (<L(2-8keV)>~10*^ erg s"^) respectively. The Seyfert 
2/Seyfert 1.8-SED class is likely to be composed by a mix- 
ture of optical Seyfert 2 and LINERs: we believe most of 
them being low luminosity AGNs, not necessarily obscured 
ones. 

In figure [2] we show an example of the IR sources whose ob- 
served SED is best-fitted by a Seyfert 2/Seyfert 1.8 template 
SED, but whose optical spectrum is typical of a spiral/star- 
forming galaxy. Note that, as discussed in detail in Section|4j 
the Seyfert 2/Seyfert 1.8 SEDs are dominated by star for- 
mation in all the MIR/FIR bands, except in the range 3-10 
^m, where the warm dust heated by the AGN shows up. 
Therefore, these objects can be considered as star-forming 
galaxies, but in our analysis we prefer to keep them as a sep- 
arate class as a reminder that they are likely to contain an 
AGN that might be important in other bands (i.e. X-rays). 

In the model we have grouped the sources into five 
"broad" SED classes, according to their SED-fitting clas- 
sification. The five populations considered in the model 
are: normal spiral galaxies (spiral: Sa, Sb, Sc, Sd, Sdm 
SEDs), galaxies powered by star-formation (starburst: 
M82, NGC6090, IRAS22491, IRAS20551, Arp220 SEDs), 
objects containing a low luminosity AGN (LLAGN: Seyfert 2, 
Seyfert 1.8 SEDs), obscured AGNs (AGN2: Markarian 231, 
IRAS19254 SEDs), and unobscured AGNs (AGNl: QSOl, 
TQSOl, BQSOl SEDs). 

For the spiral population we have assumed a SED evolv- 
ing with luminosity from an Sa spiral (Li5^m<10^ Lq) to 
an Sdm spiral (Lis^,^ ^ 10^° L©), while for the starburst 
population the SEDs vary from a moderate starburst galaxy 



(NGC6090, Li5pm==:10"' Lq), up to an extreme one (Arp220, 
Li5pm>lO^^'^L0), with fine interpolations between these 
known templates. These choices follow the conclusions of 
Gruppioni et al. (2008) that the galaxy and starburst 
SEDs evolve with 2 (and/or L), from early- to late- type and 
from moderate to extreme starbursts respectively. 
For the AGN-dominated populations (AGNl and AGN2), a 
single SED has been considered in the model as represen- 
tative of the whole class: in particular, we have assumed 
the TQSOl and Markarian 231 templates (see PoUetta et 
a. 2007 for a description) respectively. The choice of a sin- 
gle template SED as representative of each AGN-dominated 
population in the model is motivated by the fact that the 
bulk of sources (>90%) in each class is very well fitted by 
the selected templates, with apparently no clear trend with 
z or L, as instead clearly observed for spiral and starburst 
galaxies. For the same reason, for the LLAGN population we 
have assumed just two templates: the original Polletta et 
al. (2007) Seyfert 2 one, reproducing most of the ELAIS-Sl 
LLAGN, and a modified Seyfert 2 template, with a higher FIR 
bump (at A>24 ^m), needed to fit ~40% of the LLAGN SEDs 
observed by Herschel in the GOODS-N field, showing an 
excess in the FIR (see Section |3] and Gruppioni et al. 2010 
for details). The two LLAGN templates are used in the model 
in the same proportions as observed in the PEP GOODS-N 
data sample. 

As we will show in Section |4j all the template SEDs con- 
taining an AGN are "composite" SEDs, where AGN and 
starburst coexist, with the two components having different 
relative levels of importance in different wavelength ranges 
and for different SED classes. However, in the current work 
we have analysed the evolution of the different SED-classes 
as a whole, assuming that the AGN and starburst compo- 
nents influence each other and co-evolve within the same 
object. It is very hard to say how the evolution of the two 
components can differ: this would require an extremely ac- 
curate analysis and exquisite quality data to decompose the 
observed SED of each source and study the evolution of the 
AGN and starburst component separately. While such an 
analysis would be matter for a future work (Pozzi et al., in 
preparation), here we study the evolution of different SED 
populations, as resulted from a detailed SED-fitting anaysis. 
In Figure [3] we show the template SEDs used for the five IR 
populations considered in our model. 

2.2 Observational Constraints from ISO, Spitzer 
and Herschel Surveys 

The five populations considered in the model are charac- 
terised by different evolutionary properties. The evolution 
parameters are defined by minimising the differences be- 
tween data and model expectations, simultaneously for the 
single population luminosity and redshift distributions and 
for the differential source counts (total and for each pop- 
ulation) considering all the data from IR surveys available 
in the literature (from 15 ^m to 500 ^m). All the source 
counts from different surveys at the same wavelength have 
been combined by binning in flux density and averaging 
each data point weighted by its formal error (inverse of the 
squared error). The combined source counts, shown as grey 
shaded areas in Figures 5, 7, 8 and 9, have then been used 
to constrain our model. 
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The MIR data considered for constraining the global source 
counts are from several ISO and Spitzer surveys: at 15 /im 
from ELAIS-Sl (Gruppioni et al. 2002), HDF-N, HDF-S, 
Marano Firback, Marano Deep, and Marano Ultradeep (El- 
baz et al. 1999), ultradeep lensed (Metcalfe et al. 2003), 
Lockman Deep and Shallow (Rodighiero et al. 2004); at 24 
^m from GOODS (Papovich et al. 2004) and SWIRE (Shupe 
et al. 2008). The 15-fim source counts have also been con- 
strained separately for the different populations using the 
ELAIS-Sl data. The redshift and luminosity distributions 
have been constrained by using the ELAIS-Sl spectroscopic 
redshifts (La Franca et al. 2004) at 15-/im and the GOODS- 
S and -N (Rodighiero et al. 2010) and COSMOS (Le Floc'h 
et al. 2009) spectroscopic and photometric redshifts at dif- 
ferent flux density levels at 24 ^m. Using the La Franca et 
al. (2004) spectroscopic sample and the SED-based classifi- 
cation of Gruppioni et al. (2008) we have been able to sepa- 
rately constrain the redshift and luminosity distributions of 
each population at 15 fim, while at 24 fim the comparison 
has been performed on the total distributions only, since the 
24 /im samples available had no classified objects through 
SED-fitting as we would need for a direct match by popula- 
tions. 

In the FIR we have considered as constraints to source 
counts the recent results from Herschel, as well as previ- 
ous results from ISO and Spitzer. at 70 nm from Spitzer 
(CDF-S and BOOTES: Dole et al. 2004; GOODS: Frayer et 
al. 2006a; FLS: Frayer et al. 2006b; COSMOS: Frayer et al. 
2009; GOODS/FIDEL, COSMOS and SWIRE: Bethermin 
et al. 2010a) and from the shorter wavelength photometer 
(PACS; Poglitsch et al. 2010) of Herschel (PEP GOODS-S: 
Berta et al. in preparation); at 100 fj,m from 750 (Lock- 
man Hole: Rodighiero & Franceschini 2004; ELAIS: Her- 
audeau et al. 2004) and from Herschel-PACS (PEP GOODS- 
N, GOODS-S, Lockman Hole, COSMOS: Berta et al. 2010); 
at 160 /im from ISO (FIRBACK: Dole et al. 2001), Spitzer 
(CDF-S and BOOTES: Dole et al. 2004; FLS: Frayer et al. 
2006b; COSMOS: Frayer et al. 2009; GOODS/FIDEL, COS- 
MOS and SWIRE: Bethermin et al. 2010a) and Herschel- 
PACS (PEP GOODS-N, GOODS-S, Lockman Hole, COS- 
MOS: Berta et al. 2010). The global redshift distributions 
of the PEP sources at 100 /^m and 160 /tm have also been 
considered as constrain for our model. In the sub-mm we 
have considered the source counts from the recent Surveys 
with the longer wavelength instrument (SPIRE; GrifBn et 
al. 2010) on Herschel: at 250 /im, 350 /im and 500 /xm from 
HerMES (Oliver et al. 2010) and H- ATLAS (Clements et al. 
2010). 

2.3 Evolution 

For the spiral, starburst, LLAGN and obscAGN populations, 
both luminosity and density evolution arc required in order 
to fit the observables, while for unobsAGN no density evolu- 
tion is needed, just luminosity. The shape of the evolution 
functions considered here is not the commonly used (1 + 2;)'° 
with a Zbreak, but it is a more physical representation of 
the evolution, peaking at a given redshift and decreasing 
at higher redshifts, with different peaking redshifts and de- 
creasing/increasing skew rates for the different populations. 
We have modelled the evolution with z of the luminosity L 
and density p by means of a variant of the "skew-normal 



Table 2. Evolution Parameters 



SED class 


Al 


A<i, 




K 


spiral 


1.0 


2.5 


2.5 


5.0 


starburst 


12.0 


8.0 


3.5 


3.0 


LLAGN 


8.5 


6.5 


2.2 


1.8 


AGN2 


19.0 


16.0 


2.8 


0.8 


AGNl 


17.8 




4.6 


3.1 



1,4 7 




z 



Figure 4. The luminosity (solid lines) and density (daslied lines) 
evolution curves, described in Equations|2]and[3] for the five pop- 
ulations considered in the model (green: spiral, cyan: starburst, 
red: LLAGN, magenta: AGN2, blue: AGNl) that best reproduce all the 
available IR data. 

distribution" function, which is an extension of the normal 
(Gaussian) probability distribution, allowing for the pres- 
ence of skewness (i.e. Azzalini et al. 1985). In particular, we 
have considered the following shape and parametrization for 
the evolution: 

logio{L*{z)) = logio{L*{0)) + ^^e""'/""' x erf (k-) 

tJV27r V ^/ 

log,oi<f*{z)) = %io(<&*(0)) + ^l|=e-^/2-^ x erf (k-) 

Ljy Ait \ <^ / 

where erf{x) = ^ Jo '^^ "error function", Al 

and A,j> are the normalizations for the evolution of L*{z) 
and ^*{z) respectively and k (known as shape parameter, 
because it regulates the shape of the function) and u (scale 
factor) in different combinations define the skewness of the 
function and the redshift location of the peak. In order to 
limit the number of free parameters in the fit, for each popu- 
lation we have assumed the same values of k and u for both 
luminosity and density evolution. 

By considering all the available IR observational con- 
straints, we have found as best estimates for the parameters 
describing the evolution of the five populations (equations [5] 
and|3| the values reported in Table [2] The evolution func- 
tions described by these parameters and characterising the 
five IR populations of our model are shown in Figure [4] To 
allow direct comparisons, in limited redshift ranges we can 
approximate our evolution curves to those more commonly 
used, with the form (l+z)'' . We note that the spiral galax- 
ies show low evolutionary rates, both in luminosity and den- 
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Figure 5. Differential source counts normalised to Euclidean at W-fivn. Several survey data are shown (red open circles: HDF-N, green 
filled circles: HDF-S, cyan open squares: Marano Firback, asterisks: Marano Deep, magenta diagonal crosses: Marano Ultradeep, pink 
filled triangles: Lockman Shallow, orange open triangles: Lockman Deep from Elbaz et al. f999; magenta filled stars: ELAtS-Sl from 
Gruppioni et al. 2002; blue open hexagons: Abell2390 ultradeep lensed field from Metcalfe et al. 2003; blue filled triangles: Lockman 
Shallow, blue open triangles: Lockman Deep from Rodighiero et al. 2004) and compared to our model expectations. The grey shaded 
area represents the uncertainty region of the weighted averaged counts from all the different surveys. The different contributions from the 
five populations are plotted in different colours (black solid: total; green dashed: spiral galaxies; cyan dot-dashed: starburst galaxies; 
red dot-dot-dot-dashed: LLAGN; magenta dotted: AGN2; blue long-dashed: AGNl), all together {top panel), and separately by population in 
the five bottom panels (where the coloured dashed area represents the la uncertainty region of the ELAIS-Sl observed counts for that 
population) . 



sity, increasing up to 2;~0.3 (oc(l -|- z)^'^ in luminosity and 
oc(l -|- z)"'^ in density), then slowly decreasing towards the 
higher z's. Starburst galaxies evolve fast (oc(l + z)^'^ in 
luminosity and oc(l -I- z)'^'^ in density) up to z^l, with the 
luminosity and density remaining approximatively constant 
between z=l and z—2, then decreasing at higher redshifts. 
LLAGN show luminosity evolution (L(2)oc(l -I- z)'^'^) similar 
to and density evolution (p(2)oc(l + z)^'^) higher than those 
of starburst galaxies, but with a more pronounced peak at 
0~1.2-^1.4, followed by a faster decrease at z>1.5. AGNl 
evolve in luminosity at a rate oc(l -I- z)^''' up to z^l.5, then 



their luminosity remains almost constant between z~1.5 and 
z~2.5, while the evolution of the AGN2 objects shows a flat- 
tening at even higher redshift (i.e. between 2=2 and z—3), 
increasing towards the peak at a rate of «(1 -I- 2)^'^ in lu- 
minosity and ~(H- 2)^'^ in density and decreasing faster at 

2>3. 
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Figure 6. Redshift {upper panels) and luminosity {lower pan- 
els) distributions of the ELAIS-Sl 15-fim sources down to a flux 
density Si5fim=0.5 mjy. The data relative to the five popula- 
tions in the ELAIS-Sl Survey (Gruppioni et al. 2008) are shown 
as grey shaded distributions (representing the la uncertainty re- 
gion), while the overimposed model expectations are shown as 
coloured lines (same colours as in Figure pi. 



with respect to other, less statistically significant, surveys 
from smaller fields (see Figure [5|. Another little discrepancy 
between data and model consists of a high-z, high-L tail in 
the starburst modelled redshift and luminosity distribu- 
tions that is not observed in the data. We note, however, 
that we are dealing with very small numbers of objects in 
these tails: ~2 expected, versus observed. 
The 15-/^m counts are dominated by spiral galaxies and 
AGNl at SiSfim 1.5 mJy, while at fainter flux densities the 
LLAGN constitute the main population responsible for the 
peak at Sis^jm— 0.3-0.4 mJy observed in the differential 
counts. The starburst and AGN2 populations are never dom- 
inant at 15 jj,m, though their contribution is significant at 
sub-mjy level. 



3.2 The Data 

As shown in Figure [7j the agreement between data and 
model is remarkably good also at 24 /im, since the model is 
able to reproduce well the source counts, the z-distributions 
and the luminosity function in different redshift intervals 
(data from Rodighiero et al. 2010). Note that, accord- 
ing to our model, the peak observed in the counts at 
S24nm —0.3 mJy is made mainly by galaxies containing a 
LLAGN, with a minor contribution also from the AGN2 popu- 
lation. At fiux densities J 1 mJy spiral galaxies dominate 
the counts, while starburst galaxies and AGNl contribute 
almost equally and similarly to LLAGN. The redshift distri- 
bution down to S24,jm=0.15 mJy shows a peak at z~0.8-0.9 
mainly due to LLAGN and a higher- 2; tail due to the AGN2 ob- 
jects. Spiral galaxies contribute only to very low redshifts, 
while starburst galaxies peak at similar 2 to that of LLAGN, 
though the number of "pure" starburst is much lower. The 
redshift distribution shown in the plot is obtained from the 
combination of the redshift distributions in three different 
surveys (COSMOS, Le Floc'h et al. 2009 and GOODS-N 
and -S, Rodighiero et al. 2010), which were all in reasonably 
good agreement with our model. 



3 RESULTS 

3.1 The 15-^m Data 

In Figures |5] and [6] we show the results of the 15-jj.m data 
and model comparison, plotted in terms of differential source 
counts as function of flux density and of redshift and lu- 
minosity distributions (global and separately for each pop- 
ulation). Note that when comparing the single populations 
source counts, redshift and luminosity distributions from the 
model to the observed ones, we consider and apply to the 
model the same optical and 15-/im limits and spectroscopic 
incompleteness of the data. The model, constrained just up 
to the data limits, is then extrapolated to fainter optical 
magnitudes and infrared flux densities. The model is able to 
reproduce the number density, shape and peak location of 
all the distributions. A discrepancy is observed in the source 
counts, with the model overpredicting the data at flux den- 
sities of a few mJy. However, we must stress that the highly 
statistically significant SI survey, which constitutes the only 
data comparison that we used for redshifts and luminosities 
at 15 nm, in the few mJy range seems somewhat under-dense 



3.3 The 100- and 160-^im Herschel Data 

The recent results from Herschel have just provided (and 
will furtherly provide) stringent constraints to evolution 
models in the FIR, with large samples of sources with ex- 
tensive multiwavelength information, crucial not only for 
statistical studies, such as source counts, but also for de- 
tailed SED-fitting analysis and photometric redshift deriva- 
tion. The Herschel data, spanning from 70 to 500 fim, will 
systematically cover the FIR/sub-mm spectrum of thou- 
sands of galaxies, from z—Q up to 2=3-4. In particular, 
the PEP extragalactic survey samples 4 different layers with 
Herschel-PACS: from the wide and shallow COSMOS held, 
through medium size areas like Lockman Hole, to the very 
deep GOODS-N and -S areas and beyond, through grav- 
itational lensing in galaxy clusters (e.g. Abell 2218, Al- 
tieri et al. 2010). The first PEP data available, the Sci- 
ence Demonstration Phase (SDP) observations dedicated 
to the GOODS-N field, allowed us to perform a detailed 
broad-band SED-fitting analysis and to classify the 100-/im 
and 160-/xm sources using the same technique used for the 
ELAIS-Sl sources, consistently with the model populations. 
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Figure 7. Upper left panel: 24 fim differential source counts normalised to Euclidean obtained from our model (coloured curves, as 
in Fig. [sj and observed by Spitzer (magenta filled circles: Papovich et al. 2004; green filled triangles: Shupe et al. 2008; grey shaded 
area: uncertainty region of the weighted averaged counts from the two surveys.). Upper right panel: redshift distribution of the 24-/im 
sources to a flux density S24fjm=0.15 mjy: total (black line) and decomposed into populations, compared to the combined observed data 
(grey shaded area) from the COSMOS (Le Floc'h et al. 2009), GOODS-S and GOODS-N fields (Rodighiero ct al 2010). Lower panels: 
Rest-frame LP at 24 ^m in different redshift intervals: model expectations (coloured curves as in Pig. [sjl compared to observed data in 
the VVDS+GOODS fields (black circles; Rodighiero et al. 2010). 



The PEP observations provide a strong constraint to the 
SED shape of the IR sources, since for the first time the 
peak of the FIR bump has been sampled for large numbers 
of galaxies at cosmological distance. As discussed by Gruppi- 
oni et al. (2010), the template SEDs of PoUetta et al. (2007) 
generally provide very good fits to the SEDs of the PEP 
sources. However, in ~10-12% of cases the observed SEDs 
are very well reproduced by the templates over the entire 
UV/optical/NIR/MIR range, while they are systematically 
lower than the data in the FIR range. In these cases, the 
PEP 100- and 160-/im fluxes are higher by up to a factor of 
~4 than the template at the same wavelengths. This hap- 
pens mainly for the LLAGN templates (for about 40% of the 
PEP sources fitted by the Seyfert 2/Seyfert 1.8 SEDs) and 
in a smaller fraction of cases also for the spiral ones. We 
therefore constructed three new templates with a rest-frame 
0.1-15 fim spectrum similar to that of Polletta et al. (2007), 
but with a higher FIR bump, by averaging together (in 
wavelength-bins) the observed rest-frame SEDs (normalised 
to Ks band) exhibiting an excess in the FIR and fitted by 
the same template (Seyfert2/1.8/Sdm). We have therefore 



taken into account this result in our model, by considering 
both the Seyfert 2 template and the modified-Seyfert 2 one 
as representative of the LLAGN population. Since no trend 
with either luminosity or redshift has been observed in the 
PEP SDP sample for the LLAGN SEDs, in the model we have 
considered the relative contribution of the original and mod- 
ified Polletta et al. (2007) template as given by the relative 
proportion of sources fitted by the two templates. The Lumi- 
nosity Function and its evolution have then been studied for 
the five populations separately, by means of the SDP obser- 
vations, and compared to our model expectations, showing a 
very good agreement between data and model (see Gruppi- 
oni et al. 2010). The PEP SDP data in the GOODS-N have 
been supplemented with the COSMOS and Lockman Hole 
wider and shallower surveys to build galaxy number counts 
at 100 and 160 /xm (Berta et al. 2010), considered as con- 
straints to our model in the FIR. In Figure [S] we show the 
results of the comparison between the observed PEP data 
- in terms of source counts and redshift distributions - and 
our model expectations at 100 ^m and 160 ^m. 
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Figure 8. Top: Differential source counts normalised to Euclidean at 100 fim {left) and 160 fim (right). Herschel data from PEP (Berta et 
al. 2010) are represented by the red filled circles (GOODS-N), the blue filled squares (Lockman Hole) and the green triangles (COSMOS) 
in both plots. The other symbols show the previous results from ISO (100-/im: magenta filled circles from Rodighiero & Franceschini 
2004, cyan filled diamonds from Heraudeau et al. 2004; 160-fim: magenta filled circles from Dole et al. 2001) and Spitzer (160-/im: pink 
filled squares from Dole et al. 2004, orange filled triangles from Prayer et al. 2006, cyan filled diamonds from Bethermin et al. 2010a). 
The grey shaded area represents the uncertainty region of the weighted averaged counts from all the different surveys. Bottom: Redshift 
distribution of the PEP sources at 100 /im (left) and 160 /im (right) cut at a flux density corresponding to the 3cr limit (grey histogram), 
compared to model expectations to the same limit. The model curves representing the different populations arc the same as in Figures^ 
[5] and [3 



3.4 The 250-, 350- and 500-^iin Herschel Data 

In the sub-mm band the instrument BLAST performed the 
first wide and deep survey in the 250-500 nm range (De- 
vUn et al. 2009) before Herschel, producing source counts 
through a P(D) fluctuation analysys (Patanchon et al. 2009) 
and through other three different methods (e.g. 1) blind ex- 
traction using a particular algorithm called FASTPHOT; 2) 
extraction using the Spitzer 24-/im sources as a prior; 3) a 
stacking analysis of the Spitzer 24- fim galaxies in the BLAST 
data; Bethermin et al. 2010b). The HerMES (Oliver et al. 
2010) and the H- ATLAS Surveys (Bales et al. 2010) have al- 
ready observed a number of fields of different areas and sen- 
sitivity (HerMES SDP; GOODS-N, Lockman North, Lock- 
man SWIRE and FLS; H- ATLAS: 14 sq. deg. in the GAMA9 
field) using SPIRE on Herschel, providing the galaxy num- 
ber counts down to ~20 mjy at 250 fim, 350 fj,m and 500 
fim. Oliver et al. (2010) compared the observed counts with 
eight models available in the literature, showing that many 
of them cannot fit the bright end (>100 mJy) and none but 



one (Valiante et al. 2009) can fit the steep rise at 20<S<100 
mJy observed in the HerMES data. Similarly, Clements et 
al. (2010) find that none of the current models is an ideal fit 
to the H-ATLAS source counts, the best performing model 
with respect to the shape of the counts being the Negrello 
et al. (2007), but only at 500 /im, since this model overpre- 
dicts the effect of the bump at 250 and 350 fim. In Figure [9] 
(left panels) we show the comparison between our model and 
the sub-mm observed source counts. The model is able to 
reproduce very well the 250-fim and 350-fim data over the 
whole flux density range, fltting the bright end, the steep 
rise and the peak of the source counts (note the particu- 
larly good flt to the rise of the 250-/im counts, sampled with 
very high statistical significance by the H-ATLAS data). The 
HerMES local luminosity functions (LLPs) of Vaccari et al. 
(2010) are also well reproduced by our model, as shown in 
the nght panels of Figure |9] At 500 fim, we notice a slight 
discrepancy between data and model, mainly in the source 
counts, where the model tends to underestimate the data 
at 30^8500^ 100 mJy. In particular, the steep rise of the 
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Figure 9. Differential source counts normalised to Euclidean (left) and local luminosity functions (right) at the Herschel-SPIRE wave- 
lengths: 250 lira (top), 350 fim {middle) and 500 fim (bottom). The black filled circles are data from the lierMES Survey (counts: Oliver 
et al. 2010; LLPs: Vaccari et al. 2010), while the magenta filled diamonds are counts from the H-ATLAS Survey (Clements et al. 2010). 
The open symbols are source counts derived from BLAST data by Bethermin et al. (2010b; orange squares: deep field with and without 
24-/im priors; blue circles: counts computed with a stacking analysis on 24-fim positions) and Patanchon et al. (2009; red stars; P(D) 
analysis). The grey shaded area represents the uncertainty region of the weighted averaged counts from all the different surveys. 



source counts starts at brighter flux densities (~100 mjy) in 
the data than in the model, the latter predicting a smoother 
rise up to ~20-30 mJy. The model LLF at 500 /xm is flat- 
ter than the data in the two lower luminosity bins, though 
the knee and the bright tail of the local luminosity func- 
tion is well reproduced. The discrepancy between data and 
model observed at 500 fim could be due to the presence 
(even locally) of a "cold" population, whose SEDs are not 
well represented by the templates considered by our model 



(which are just extrapolated in the sub-mm regime, not us- 
ing the SPIRE and BLAST data as constraints). Indeed, 
z<l analogues of the "cold" z>l sub-mm galaxies have al- 
ready been found by Symeonidis et al. (2009) before Her- 
schel, while Elbaz et al. (2010) find lower dust temperatures 
than previously inferred (due to the lack of constraints at 
FIR wavelengths before Herschel) in z<l-1.5 galaxies de- 
tected by PA CS' and SPIRE. Moreover, they also find SEDs 
"colder" than those of their local analogues for a significant 
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fraction (10-20%) of LIRGs and 2~1.6 ULIRGs. On the 
other hand, Shultz et al. (2010), by analysing the colours of 
the HerMES sources, discovered a population of red bright 
objects that may consist mostly of "colder" SEDs, but with 
a fraction (>12%) of distant lensed ones. We would there- 
fore need to further investigate the FIR/sub-mm SEDs of 
the IR populations: as soon as we could have access to the 
SPIRE data for large samples of galaxies, we will be able to 
properly fit the observed SEDs up to 500 ^m, comparing, 
adjusting and modifying our templates to better reproduce 
the real Universe. In addition, since a considerable fraction 
of sub-mm bright sources are expected to be lensed by fore- 
ground galaxies (e.g. according to Negrello et al. 2007 all the 
500-/xm sources brighter than 100 mjy and with 2<z<3 are 
lensed; see also the recent H-ATLAS results of Negrello et 
al. 2010), the effect of lensing should be properly taken into 
account when performing statistical studies (with both data 
and models) like source counts and luminosity functions. 



4 THE AGN CONTRIBUTION TO THE 
INFRARED EMISSION 

As shown in the previous sections, from our model we ex- 
pect a significant contribution to the source counts and lu- 
minosity functions from objects containing an AGN. How- 
ever, AGNs and starbursts often co-exist in the same object, 
while in the IR populations defined by our different SEDs, 
the real AGN contribution is not disentangled from that due 
to star-formation. Here we try to identify - in a very sim- 
plified way - the AGN contribution to the monochromatic 
and total IR luminosity emitted by the different populations 
considered in our model. To this purpose, we have decom- 
posed, by mean of a standard minimisation technique, 
the template SEDs into three distinct components: a stellar 
component emitting most of its power in the optical/NIR, 
an AGN component due to hot dust heated by the central 
black hole, and peaking in the MIR, and a starburst (SB) 
component representing the major contribution to the FIR 
spectrum. The algorithm combines synthesis stellar mod- 
els built using the Padova evolutionary tracks (Bertelli et 
al. 1994), AGN dusty tori models from Fritz et al. (2006) 
and 6 empirical starburst SEDs to reproduce the observed 
broad-band spectra. We refer to Pozzi et al. 2010 (see also 
Hatziminaoglou et al. 2008, 2010; Vignah et al. 2009) for 
a detailed description of the method and of the individ- 
ual model components. In particular, we run the algorithm 
on the template SEDs representative of our AGN popula- 
tions (LLAGN, AGN2and AGNl): Seyfert 2, Markarian 231 and 
TQSOl of PoUetta et al. (2007), respectively 

In Figure [lO] we show the results of the template de- 
composition, with the AGN contribution highlighted in red. 
Note that the optical/NIR emission in the LLAGN spectrum 
is almost exclusively due to the stellar component from the 
host galaxy, while in the AGNl spectrum the emission from 
the accretion disk is dominant in that wavelength range. The 
AGN2 SED shows an intermediate situation between LLAGN 
and AGNl, though apparently closer to that of the AGNl, with 
the torus emission dominating in the whole MIR range. How- 
ever, we note that the torus emission in the AGN2 spectrum, 
at odds with the AGNl one, must be of "type 2" (i.e. edge- 
on), such as to reproduce the strong silicate absorption at 
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Figure 10. Templates representative of the three populations 
containing an AGN {left: LLAGN; middle: AGN2; right: AGNl) de- 
composed into three different contributions; a stellar component 
emitting most of its power in the optical/NIR (green), an AGN 
component due to hot dust heated by the central black hole, and 
peaking in the MIR (red), and a starburst representing the ma- 
jor contribution to the FIR spectrum (blue). The total SED is 
represented by the black solid fine. 

9.7 /^m (typical of heavily obscured sources) observed in the 
IRS spectrum of our prototypical AGN2 object Markarian 
231 (Weedman et al. 2005). Though Markarian 231 is a well 
known heavily obscured BAL QSO, the Mrk231 (AGN2) tem- 
plate reproduces the observed SED of "obscured" AGNs re- 
gardless of their optical spectra (i.e. both type 1 and type 2 
in the optical; Gruppioni et al. 2008) and also those of many 
Dust Obscured Galaxies (DOGs; Dey et al. 2008; Lanzuisi 
et al. 2009), with extreme MIR-to-optical ratios (Gruppi- 
oni, Vignali, Fritz et al. in preparation), likely to harbour 
obscured AGNs. The modelled 9.7-/xm feature is in absorp- 
tion in the LLAGN SED, while it is in emission in the AGNl 
one, as commonly observed in type 2 and 1 AGN spectra 
respectively (but see i.e. Mason et al. 2009 for examples 
of silicate emission in type 2 AGNs). Using these SED de- 
compositions, we have derived the fractional contributions 
of AGN, starburst (SB) and evolved stars to the monochro- 
matic and bolometric (both 8-1000 and 1-1000 ^.m) 
IR luminosities for the three representative template SEDs. 
These fractions, reported in Table [S] have been used to dis- 
entangle the contribution to the total IR luminosity function 
due to SB and AGN activity. The evolved stars component in 
the decomposition of the AGN2 SED was necessary to repro- 
duce the optical/NIR part of the spectrum, since the "type 
2" torus needed for modelling the silicate feature in absorp- 
tion, in the Fritz et al. (2006) model does not account for 
emission observed in these bands. Different tori models (i.e. 
the "clumpy" one of Nenkova et al. 2008) can explain the 
optical spectra of type 2 AGNs in terms of scattering con- 
tribution with no need for a stellar contribution. However, 
whether there is or not a stellar component in the AGN2 tem- 
plate is not relevant for our further analysis and conclusions 
and is beyond the scope of this work. 

4.1 Star- Formation and Accretion History 

We have used the LFs provided by our model, and the SED 
decomposition described above, to estimate the cosmic evo- 
lution of the star-formation density (SFD), psfr{z) (de- 
rived from the total comoving IR luminosity density due to 
star-formation, P^r,{z)), and of the SMBH accretion density, 
^ bhar{z), as a function of redshift. The comparison be- 
tween these two quantities is a crucial tool for understanding 
galaxy and AGN evolution and the role played by the AGNs 
in the formation of galaxies. We have estimated Pjr{z) as 
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Table 3. AGN, SB and stellar contribution to the IR luminosity 
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Figure 11. Redshift evolution of the total IR luminosity density 
pj^{z) (and SFD, Psfr) ^ expected from our model (yellow 
shaded area) and compared with estimates from IR surveys (green 
shaded area: Le Floc'h et al. 2005; orange filled triangles: Caputi 
et al. 2007; blue filled squares: Rodighiero et al. 2010; black filled 
circles: Gruppioni et al. 2010, with the arrows standing for lower 
limits due to luminosity incompleteness). The contributions to 
PiR (^'^ Psfr) from the different populations are shown as green 
dashed line, cyan dot-dashed line, red dot-dot-dot-dashed lines, 
magenta dotted line and blue long-dashed line for the spiral, 
starburst, LLAGN, AGN2 and AGNl class respectively. We have cal- 
culated the upper and lower limits of our total SFD expectation 
by considering as AGNs only those spectroscopically confirmed 
and all the SED-fitting derived ones respectively (see Gruppioni 
et al. 2008). On the right y-axis we report the SFD obtained by 
means of the Kennicutt (1998) conversion and, for comparison, 
we plot (as grey shaded area) the evolution of the SFR density 
Psfr obtained from a large compilation of extinction-corrected 
UV/optical and Ha data by Hopkins & Beacom (2006). Given 
the lack of data constraints at high-z, our model expectations are 
highly uncertain at z>2.5. 



follows: 



SF , ^ 
PiR (Z) 



o4>{Ls 



))dlogLs 



(4) 



where Lfi'iooo is the 8-1000 fim IR luminosity due to star- 
formation. The SFD has been derived from the IR lumi- 
nosity, according to the conversion of Kennicutt (1998): 
PSFfl=1.7xlO"^° piR [Mq yr"^ Mpc"^] (for a Salpeter 
IMF). 

we show Piji{z) 



In Figure 
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as expected from our 
model and compared with estimates from different IR sur- 
veys (i.e. at 24 ^m: Le Floc'h et al. 2005; Caputi et al. 
2007; Rodighiero et al. 2010; at 100 p,m and 160 pm: Grup- 
pioni et al. 2010). To derive the uncertainties of our cal- 
culations, we have considered as upper limit for the AGNs 
the number of AGNs derived from the SED fitting analysis 
and as lower limit the fraction of AGNs spectroscopically 
confirmed through optical emission lines (see La Franca et 
al. 2004; Gruppioni et al. 2008). By considering these val- 
ues we obtain the lower and the upper limit to the esti- 
mate of pf^{z) respectively. We have also compared the 
SFD evolution obtained from the IR to that derived from 
optical/UV and Ha observations, presented in a large and 
homogenised collection from different surveys by Hopkins 
& Beacom (2006). Our model expectation is in very good 
agreement with all the IR data estimates, confirming the 
rapid increase of pf^ up to z^l. The increase of the IR lu- 
minosity density is followed by a peak at l<z<2 and by a 
decrease from 2~2-2.5 towards the higher redshifts. At the 
lower redshifts (2:<0.3) the IR luminosity density is domi- 
nated by the spiral population, while in the 0.3<2<2-2.5 
range, the principal contributors to pf^ are galaxies with a 
LLAGN SED. Pure starburst galaxies are also important in 
the same redshift interval, but are never dominant at any z: 
at z>2.5 the AGN2 SED objects start dominating up to 2>4, 
when they are overtaken by the AGNl population. Therefore, 
galaxies likely to contain a low-luminosity AGN and galax- 
ies powered by a starburst are responsible for the peak of 
the IR luminosity density at 2~l-2, then galaxies hosting 
increasingly powerful AGNs become increasingly important 
towards the higher z's. 

The total IR emissivity predicted by our model at z>2 
is lower and shows a faster convergence than previously pub- 
lished results based on UV/optical or Ha observations (i.e. 
Hopkins & Beacom, 2006), which are subject to large ex- 
tinction corrections, or on 24-/im data (i.e. Perez-Gonzalez 
2005), which require large spectral extrapolations to derive 
the total IR luminosity, especially at high redshifts. Sim- 
ilar high-z convergence of the total IR luminosity density 
has been found recently by Franceschini et al. (2010) with 
a model based on the analysis of a large IR database of 
high-redshift galaxies at long wavelengths. However, up to 
now the total emissivity of IR galaxies at high redshifts is 
poorly constrained, due to the scarcity of Spiteer galaxies at 
z>2 and the incomplete information on the ^-distribution 
of sub-mm sources (Chapman et al. 2005). At high redshifts 
also our model expectations are very uncertain, since no con- 
straint from data at 2:>2-2.5 are available from pre-Herschel 
Surveys. The large numbers of high-z galaxies provided by 
the deep Surveys performed with Herschel (and furtherly 
with SPICA) will be crucial to assess galaxy and AGN evo- 
lution in the IR at z>2.5. 
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Figure 12. Redshift evolution of the BHAR density, bhar{z), 
estimated from our IR analysis. Pink shaded area (solid lines) 
shows the total accretion rate, while the red (dot-dot-dot-dashed 
lines), blue (long dashed lines) and magenta (dotted lines) show 
the accretion rate due to LLAGN, AGNl and AGN2 objects respec- 
tively. The uncertainty regions for each population have been ob- 
tained by considering AGNs all the SED-fitting derived ones and 
only those spectroscopically confirmed (see Gruppioni et al. 2008) 
in the calculation of the upper and lower envelopes respectively. 



The SMBH accretion rate can be derived once the 
bolometric luminosity function z) is known, where 

^boi^ ~ ^radMc? is the iutriusic bolometric luminosity pro- 
duced by a SMBH accreting at a rate of M with a radia- 
tive efficiency trad- Crucial factors in the determination of 
$(L;f„f^,2) (Hopkins, Richards & Hernquist 2007; Merloni 
& Heinz 2008) are the completeness of any AGN survey 
and a suitable correction to estimate, from observations in 
one band, the AGN bolometric luminosity (i.e. bolometric 
correction BC). Observations in the hard X-ray band are 
commonly used, since the X-ray surveys provide a relatively 
unbiased census of AGN in the Universe. Here we estimate, 
for the first time, the SMBH accretion from the IR emission 
originating from the circumnuclear dusty material that in- 
tercepts a fraction of the inner optical-UV thermal accretion 
disk emission. While the nature of the dusty material is still 
a matter of debate (i.e. smooth vs. clumpy distribution), 
dust around the BH is observed in almost all the AGNs 
(see the review of Elitzur 2008 and references therein). As 
bolometric correction BC we used a value of ~1.5 for the 
AGNl and AGN2 templates, and BC~2 for the LLAGN template. 
These bolometric corrections are first-order empirical esti- 
mates derived in Pozzi et al. (2007), where the broad- band 
SEDs of a sample of X-ray selected AGN have been stud- 
ied. As explained in Pozzi et al. (2007), these corrections 
take into account the geometry of the torus (i.e. covering 
factor /) and the effects of orientation (i.e. effective opti- 
cal depth rg.T^m of the dust along the line of sight). The 
geometry correction is based on statistical arguments, em- 
ploying the ratio between obscured and unobscured quasars 
as required by recent X-ray background synthesis models 
(Gilli et al. 2007) to infer a typical torus covering factor of 
/ 0.67 (hence i?C=1.5). As reported by Vasudevan et al. 
(2010), this value is also consistent with the covering fac- 
tor obtained from recent detailed clumpy torus models (i.e. 
Nenkova et al. 2008) for typical torus parameters (e.g. num- 
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Figure 13. Redshift evolution of the BHAR density ^bhar{^)^ 
multiplied by a constant factor (equal to 500) is plotted as a pink 
shaded area alongside to our SFD derivation shown in Figure [TT] 
(yellow shaded area) and to the best fit to a large compilation 
of optical/UV data for the SFD (grey shaded area, from Hopkins 
and Beacom 2006). 

ber of line-of-sight clouds ~ 5, opening angle ~30-45°). The 
effect of orientation was computed from the ratio between 
the luminosities of a face-on and an edge-on AGN using the 
AGN SEDs of Silva et al. (2004). Consistent results were 
derived in Pozzi et al.(2010) using a more sophisticated ap- 
proach, i.e. computing the bolometric corrections using the 
best-fitting torus models, as the ratio between the accre- 
tion disk luminosity given as input to the radiative transfer 
model and the observed reprocessed infrared luminosity in 
output (see Fig. 4 from the pioneering work of Pier & Krolik 
1992). Since the latter approach implies an exploitation of 
the degeneracy of the torus solutions which is beyond the 
scope of the present work, we prefer the first, simpler and 
straightforward, method. 

The expression used for estimating the SMBH accretion 
(BHAR) is: 

,T, t \ [ — ^rad) BC J,/ r \ J7 r 

9bhar(z) = / ^ (t>(Li^iooo)dlogLi^ 

Jo ^radC 

where BC is the bolometric correction to the 1-1000 /im 
IR luminosity depending on the SED type and L^^-^qq indi- 
cates the 1-1000 nm IR luminosity due to the AGN. For erad 
we have assumed the canonical value of 0.1 (see e.g. Hop- 
kins, Richards & Hernquist 2007), but changing this value 
would simply result in a change of the overall normaliza- 
tion. In Figure [12] we show the predicted 'i^BHARiz), with 
the contribution from the different AGN populations. The 
uncertainty region of the single contributions and of the total 
estimate have been obtained by considering, as for the SFD, 
the SED-fitting and spectroscopic fractions of the different 
AGN types respectively as upper and lower limits in our cal- 
culations. We note that for the AGNl and AGN2 SED objects 
the difference between the numbers derived by SED-fitting 
and by optical spectroscopy is small, while for the LLAGN the 
fraction of spectroscopically confirmed AGNs is only about 
30% of the SED-fitting one (see Gruppioni et al. 2008). How- 
ever, even considering the lower limits for objects containing 
an AGN, which conservatively takes into account the possi- 
bility that a significant fraction of our LLAGN are not AGNs, 
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but starburst galaxies not well reproduced by our limited 
set of templates, the global result of our analysis does not 
change significantly. 

The AGNl contribution dominates the BHAR, especially at 
low (2:<0.3) and high redshifts {z>3), with the LLAGN and 
the AGN2 population's accretion density peaking at l<z<2 
and 2<z<3 respectively and reaching the AGNl accretion val- 
ues just in these redshift ranges. The BHAR obtained from 
our IR estimate is reasonably consistent in shape with previ- 
ous derivations from X-rays (i.e. Merloni, Rudnick & di Mat- 
teo 2004; Merloni & Heinz 2008), though it is about a factor 
of ~2 higher at 1 < z 2 3, where the BHAR peaks. Note that 
the derivation of BHAR from IR luminosity is completely 
independent from that obtained from X-ray data and re- 
quires bolometric corrections at least a factor of ~10 lower 
and probably less uncertain than those required from obser- 
vations in the X-ray band (~1.5-2 against ~10-40 in the 
X-rays). Moreover, to derive the BHAR from X-rays, sub- 
stantial assumptions need to be made regarding the number 
and redshift distribution of compton-thick AGNs missed by 
the X-ray surveys. These obscured objects should instead be 
observable in the IR and are, in principle, already included 
in our calculations. The main source of uncertainty in our 
approach is due to the identification of sources containing a 
LLAGN: in our calculations we have considered as lower value 
that obtained by considering as AGN only those spectroscop- 
ically confirmed. Only future space missions with high reso- 
lution spectrometers in the MIR and FIR range, like MIRI 
(Wright et al. 2004) on JWST and SAFARI (Swinyard et 
al. 2009) on SPICA, will be able to definitely confirm or 
reject the presence of low-luminosity AGNs inside objects 
classified as LLAGN on the basis of their SEDs. 

The value of local BH accreted mass that we 
obtain by integrating our 'i' bhar{z) over time (i.e. 
PBH,o=J^ •^BHAR{z){dt'/dz)dz, with dt/dz^l/[Ho{l + 

z)\/f^m(l + z)^ + Qa] being the differential cosmic time as 
a function of redshift) is psH,o=6.5-9.2xlO'^ M© Mpc'^ 
Despite the large uncertainties and the simplifying as- 
sumptions in our calculation, we find that our value of 
pBHfi is in broad agreement with (although somewhat 
on the high side of) previous derivations from X-rays, 
of QxlO'^ M0 Mpc"^ {Q.l/trad) (Barger et al. 2001), 
of (7.5-16.8) X 10^ (0.1/e^ad) M© Mpc"^ (Elvis, Risaliti & 
Zamorani 2002), 4.811^ 99 x 10=^(0. l/e,.<,<i) M© Mpc"^ (Hop- 
kins, Richards & Hernquist 2007) and (3.2-5.4) x 10^ Mq 
Mpc~^ (Shankar, Weinberg & Miralda-Escude 2009). 



5 PREDICTIONS FOR FUTURE SURVEYS 
WITH JWST- MIRI AND SPICA-SAFARI 

We have used the model described above to make predictions 
for surveys to be performed with the MIRI and SAFARI in- 
struments on board of the JWST and of the JAXA-ESA 
satellite SPICA respectively. For our simulations, we chose 
two representative bands of MIRI (i.e. 10 and 18 pm) and 
two of SAFARI (i.e. 
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we show 



and 85 ^m). In Figure 
the extragalactic differential source counts and the redshift 
distributions expected in the four bands. The adopted flux 
limits in our simulations are the estimated lOcr detection 
level for a lO** seconds exposure with MIRI at 10 and 18 fim 
(0.7 /iJy and 4.3 /iJy respectively, as provided by the MIRI 



science team at the official MIRI web pages) and close to the 
estimated 5^ _R4 7?/ confusion limits at 48 and 85 /xm (~0.015 
mjy and ~0.5 mjy respectively, estimated by means of our 
model and based on the telescope specification given by the 
SPICA Study Team Collaboration, 2009, and on the recent 
proposed rescope of the SPICA telescope, as the change of 
the primary mirror diameter from 3.5 to 3.2 m). To these 
limits, we expect to detect ~6xl0* sources/deg^ at 10 /xm 
and ~1.5xl0^ sources/deg^ at 18 with MIRI, while 
~2xl0^ and -8x10" sources/deg^ with SAFARI &t 48 and 
85 /im respectively. 

The deep Surveys that will be performed with SAFARI at 
85 /im (i.e. 60-110 /im band) and 160 p,ra (i.e. 110-210 yLva 
band) will easily reach the confusion limits. In particular, 
given the nominal field of view of SAFARI of 2' x 2', at 
85 /im it will be possible to observe down to confusion a 
field as large as one of the GOODS areas (~ 10' x 15') in 
just 2.5 minutes and a COSMOS-like field (2 deg^) in ~2 
hours (not considering the overheads). On the other hand, 
at 48 /im (i.e. 35-60 /im band) about 26 hours would be 
needed to cover one field like GOODS and ^ 1000 hours 
for a larger field like COSMOS to 0.015 mJy (timely pro- 
hibitive) . These confusion limited surveys with SAFARI will 
be > 3—10 times deeper than the deepest GTO Survey with 
Herschel PACS (i.e. PEP 100 /im Survey in the GOODS-S 
to 1.7 mJy), detecting, in a GOODS-size field, hundreds of 
L7ij<10^'' L0 galaxies at z^l and few hundreds of LIGs at 
z~2. In addition, a confusion limited survey in the 35-60 
/im band will detect also of the order of 150 AGNs (all with 
i/i? > 10" Lq) at 2>3. The 110-210 /im will get confusion 
limited very soon (estimated Scon/(5(7)~13 mJy), therefore 
this band is more suited for large and moderately shallow 
surveys to detect rare and luminous objects. 

From Fig. [14] we can notice how the dominating popu- 
lations are expected to change with increasing wavelength, 
with the redshift distribution at 10 /tm being largely domi- 
nated by spiral galaxies and that at 85 /im being dominated 
by Seyfert 2-SED objects, with a major contribution of spi- 
ral galaxies at lower redshifts (z;^0.5). At 10 /im we expect a 
significant fraction of spiral galaxies also at high-z (between 
z^2 and 4), not present in the other bands. 

By considering the limiting fiuxes discussed above, we 
have computed the limiting luminosities for the prospected 
Surveys with MIRI and SAFARI and estimated the evolu- 
tion of the SFD and BHAR that we expect to obtain with 
real data from these surveys (using the same procedure de- 
scribed in Section [4.1 1 but integrating the total IR LF only 
down to the limiting luminosities). In Figure 15 we show 



the SFD and BHAR (multiplied by a factor of 500) that we 
expect to obtain with data from the confusion limited Sur- 
vey with SAFARI at 48 /tm, compared to the predicted SFD 
and BHAR from our model, corresponding to an "ideal" FIR 
survey covering all the luminosities at all redshifts (the same 
curves shown in Fig. 13 1. For comparison, in the upper panel 



of Figure |15| we show the same simulation with the deepest 
PEP Survey in the GOODS-S, reaching 1.7 mJy at 100 /tm 
(the GOODS-Herschel Survey will reach ~0.6 mJy, but just 
in a very small area of 50 arcmin^). In agreement with the 
results of Gruppioni et al. (2010; see also Fig. Ill in the 
GOODS-N, our predictions show that PEP is complete in 
SFD up to 2:~1-1.5, becoming more and more incomplete 
with increasing redshift (i.e. at 2=3 the incompleteness in 
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Figure 14. Differential extragalactic source counts normalized to tlie Euclidean slope (left) and redshift distributions {right) as predicted 
by our model in four selected bands that will be covered by future Surveys with MIRI (i.e. 10 and 18 fim) and SAFARI (i.e. 48 and 85 
/xva). The redshift distributions have been simulated to the estimated point source lOo" detection level for a 10"* seconds exposure with 
MIRI aX 10 and 18 ^m (provided by the Mffl/ science team at the official MIRI -weh pages) and close to the estimated 5y4R4_R/ confusion 
limits (given the specification of the SPICA Study Team Collaboration, 2009, and the recent proposed rescope of the telescope) at 48 
and 85 fim. These limits are shown by vertical dashed lines in the left panels. 
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Figure 15. BHAR density, '^bhar (blue vertical-dashed area), 
and SFD, Psfr (green diagonal-dashed area), as predicted by our 
model for the deepest PEP Herschel Survey at 100 fim (top) and 
for a confusion limited Survey with SAFARI at 48 fira (bottom) 
and compared to the "total" model expectations reported in 
Fig. [13] 

SFD due to the survey fiux limit is about a factor of 3 or 
more) . The BHAR that we could measure with PEP is com- 
plete up to z~1.5-2, then decreasing less rapidly than the 
SFD up to z—3 and dropping down at higher ^'s. With the 
SAFARI Survey we expect to be able to measure almost all 
the SFD to and most of it to 2~3-3.5, and almost all 
the BHAR to 2;~3. Moreover, as mentioned in the previous 
Section, the high resolution spectrometer of SAFARI will be 
crucial in identifying AGNs and separating the SB from the 
AGN contribution, to measure with great precision the SFD 
and BHAR in the high redshift Universe. 



6 CONCLUSIONS 

We have presented a new model for galaxy and AGN evo- 
lution in the IR, in which the evolutionary properties of 
the different IR populations, defined by means of a detailed 
SED-fltting analysis (see Gruppioni et al. 2008; 2010), are 
separately studied and constrained by all the available re- 
sults from MIR and FIR surveys. The model identifies five 
main SED classes: three containing an AGN at different lev- 
els of dominance (AGNl where the unobscured AGN dom- 
inates up to the MIR domain, LLAGN with a galaxy SED 
but showing a flattening at MIR wavelengths explainable 



only by the presence of a dusty torus, and AGN2 where an 
obscured AGN and starburst co-exhist) and two powered 
by star-formation only (normal spiral and moderate-to- 
extreme starburst galaxies). 

The main results of this work are summarised as follows: 

• The spiral galaxies show low evolutionary rates, both 
in luminosity and density, increasing up to 2:~0.3 (oc(l-|-2)^'^ 
in luminosity and oc(l + z)"'* in density), then slowly de- 
creasing towards the higher z's. Starburst galaxies evolve 
fast (oc(l + z)"^'^ in luminosity and oc(l -I- z)^''' in den- 
sity) up to z'^1, with the luminosity and density remain- 
ing approximatively constant between 2=1 and z=2, then 
decreasing at higher redshifts. LLAGN show luminosity evo- 
lution (L(z)cx(l + z)^'^) similar to, and density evolution 
{p{z)cc{l + z)^'^) higher than, those of starburst galaxies, 
but with a more pronounced peak at 2:~1.2-^1.4, followed 
by a faster decrease at z>1.5. The AGNl luminosity evolves 
as (x(l -I- z)^'^ up to z~1.5, then remains almost constant 
between z~1.5 and 2~2.5, while the evolution of the AGN2 
objects shows a flattening at even higher redshift (i.e. be- 
tween z—2 and z=3), increasing towards the peak at a rate 
of ^(l + z)'^'"' in luminosity and ^(l + z)'^'^ in density and de- 
creasing faster at z>3. With our backward evolution model 
we are able to reproduce well the source counts, redshift 
distributions and luminosity functions in the MIR and FIR 
bands, provided by all the main IR space observatories (i.e. 
ISO, Spitzer and Herschel). 

• We have decomposed the template SEDs representative 
of the populations containing an AGN into three distinct 
components: a stellar component emitting most of its power 
in the optical/NIR, an AGN component due to hot dust 
heated by the central black hole, and peaking in the MIR, 
and a SB component representing the major contribution 
to the FIR spectrum. In this way, we have estimated - al- 
though in a very simplifled way - the AGN contribution to 
the monochromatic and total IR luminosity emitted by the 
different populations considered in our model. 

• Using the LFs given by our model and the total IR 
luminosity due to SF derived from our template SED de- 
composition, we have estimated the cosmic evolution of the 
total IR luminosity density due to star- formation, pfjiiz), 
as a function of redshift. Our model expectation is in very 
good agreement with all the IR data estimates, confirming 
the rapid increase of pf^ up to 2~1. The IR luminosity den- 
sity shows a peak at l<z<2 and a decrease from z~2-2.5 
towards the higher redshifts. At z<Q.3 pj^ is dominated 
by the spiral population, while in the 0.3^2^2-2.5 range 
the principal contributors are galaxies with a LLAGN SED. 
Starburst galaxies are also important in the same redshift 
interval, but are never dominant at any z. The AGN2 SED 
objects start dominating at z>2.5 up to 2~4, when they 
are overtaken by the AGNl population. Star-forming galaxies 
containing or not a low-luminosity AGN (with LIRG lumi- 
nosities) are therefore responsible for the peak of the IR 
luminosity density at z~l-2, then galaxies hosting increas- 
ingly powerful AGNs (in the ULIRG luminosity range) be- 
come increasingly important towards the higher z's. 

• For the first time, the SMBH Accretion Density 
'i' bhar{z) as a function of redshift has been derived from 
IR rather than X-ray data. The BHAR obtained from our 
IR estimate is reasonably consistent in shape with previous 
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derivations from X-rays (i.e. Merloni, Rudnick & Di Mattco, 
2004; Merloni & Heinz 2008), though it is slightly higher at 
liziS, where the BHAR peaks. The Compton-thick AGNs 
are already included in our calculations, while substantial 
assumptions need to be made regarding their number and 
redshift distribution in the X-rays. The AGNl contribution 
dominates the BHAR, especially at low («<0.3) and high 
redshifts {z>3), with the LLAGN and the AGN2 population's 
accretion density peaking at l<z<2 and 2<z<3 respectively 
and reaching the AGNl accretion values just in these redshift 
ranges. 

• We have simulated source counts, redshift distributions 

and SFD and BHAR that wc expect to obtain with the 
future cosmological Surveys in the MIR/FIR that will be 
performed with JWST-MIRI and SPICA-SAFARI. 
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